A Reinvestigation of Moving Punctured Black Holes with a New Code 



Zhoujian Cao* 

Institute of Applied Mathematics, Academy of Mathematics and System Science, 
Chinese Academy of Sciences, Beijing 100080, China, 
Theoretical Institute for Advanced Research in Astrophysics, Academia Sinica, Taiwan, 
State Key Laboratory of Scientific and Engineering Computing, China 

Hwei-Jang Yo^ and Jui-Ping Yu 
Department of Physics, National Cheng-Kung University, Tainan 701, Taiwan 
(Dated: December 4, 2008) 

We report on our code, in which the moving puncture method is applied and an adaptive/fixed 
mesh refinement is implemented, and on its preliminary performance on black hole simulations. 
Based on the BSSN formulation, up-to-date gauge conditions and the modifications of the formu- 
lation are also implemented and tested. In this work we present our primary results about the 
simulation of a single static black hole, of a moving single black hole, and of the head-on collision of 
a binary black hole system. For the static punctured black hole simulations, different modifications 
of the BSSN formulation are applied. It is demonstrated that both the currently used sets of mod- 
ifications lead to a stable evolution. For cases of a moving punctured black hole with or without 
spin, we search for viable gauge conditions and study the effect of spin on the black hole evolution. 
Our results confirm previous results obtained by other research groups. In addition, we find a new 
gauge condition, which has not yet been adopted by any other researchers, which can also give stable 
and accurate black hole evolution calculations. We examine the performance of the code for the 
head-on collision of a binary black hole system, and the agreement of the gravitational waveform it 
produces with that obtained in other works. In order to understand qualitatively the influence of 
matter on the binary black hole collisions, we also investigate the same head-on collision scenarios 
but perturbed by a scalar field. The numerical simulations performed with this code not only give 
stable and accurate results that are consistent with the works by other numerical relativity groups, 
but also lead to the discovery of a new viable gauge condition, as well as clarify some ambiguities in 
the modification of the BSSN formulation. These results demonstrate that this code is reliable and 
ready to be used in the study of more realistic astrophysical scenarios and of numerical relativity. 

PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf, 97.60. Lf 



I. INTRODUCTION 

There are two main purposes for studying numerical 
relativity: One is to investigate the mathematical issues 
of the geometry of the Einstein manifold. These issues 
include cosmic censorship, the hoop conjecture, the Pen- 
rose inequality, and critical phenomenon The other, 
more practical one, is to study the dynamics of astro- 
physically compact objects. To meet the needs of existing 
(e.g., LIGO [2|, VIRGO H, GEO600 [|, and TAMA [|) 
and planned (e.g. LISA Q) gravitational wave detectors, 
the theoretical prediction of the gravitational waveform 
for realistic sources has become urgent. In both of the as- 
pects of numerical relativity, the most important targets 
for study are black holes and neutron stars, especially 
those in binary systems. 

After decades of exploration by many researchers, re- 
cently there have been exciting breakthroughs in the sim- 
ulation of the evolution of binary black holes (BBHs) 
0i B @- Soon after these breakthroughs, the moving 
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puncture method based on the BSSN formalism [KJl was 
widely used by numerical relativity groups to deal with 
black hole systems [ll|, [H, and neutron star systems 
[l^ . Much interesting physics related to BBHs has been 
explored, including the waveform of the gravitational ra- 
diation [3, Ell, the spin-orbit coupling effect and 
the recoil velocity [iSl [Toj . It has been emphasized in 
almost all of these studies that the treatment of the sin- 
gularity problem, i.e., the moving puncture method and 
the gauge condition, are both key factors in this series 
of successes. Not needing any inner boundary condition 
makes the moving puncture method much simpler to han- 
dle, compared with the excision method [20] . This sim- 
plicity has made the moving puncture method popular 
in the numerical relativity community. The Bona-Masso 
type slicing gauge conditions [21] for the lapse function 
and many driver gauge conditions (e.g., the F-driver) for 
the shift vector [22, [23| have been shown to be very im- 
portant in making the moving puncture method work. 
However, in it was shown that, although the details 
of the gauge conditions used in the punctured BBH evo- 
lutions are different, only certain gauge choices allow one 
to evolve a single moving puncture black hole. It is desir- 
able to better understand the effect of the gauge choices 
on black hole evolutions. 
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In order to study numerical relativity and also investi- 
gate the specific aforementioned topics we develop, from 
scratch, a code based on the moving puncture method. 
We adopt a fourth-order finite differencing scheme for 
the spatial derivatives and the Crank-Nicholson scheme 
for the time integrations. We use the GrACE package 
[23 | to implement both the mesh refinement and the par- 
allelization in our code. In this paper we present our pri- 
mary results about the evolution of a single static black 
hole, of a single moving black hole, of the head-on col- 
lision of a BBH, and of the head-on collision perturbed 
by a scalar field. Our results confirm many results ob- 
tained in the previous works by other groups. On the 
other hand we find a new gauge condition, which has not 
been tried by other researchers, that can also give stable 
and accurate black hole evolution calculations. We also 
observe the effects of the existence of a massless scalar 
field in delaying the head-on collision, depending on the 
initial configuration of the scalar field. All of these results 
enhance our confidence in this code, and thus we will ap- 
ply the code to more realistic astrophysical calculations 
in the near future. 

The remainder of the paper is organized as follows: 
In Section ini we first summarize the BSSN formulation, 
the conventional modifications and adjustments, and the 
gauge conditions. Then we describe the numerical meth- 
ods used in this code in section Hill including the details 
of the FMR/AMR algorithm, the finite-differencing sten- 
cils, and the Kreiss-Oliger dissipation. In Section Hvl the 
initial data is outlined. In Section |Vl we present our nu- 
merical results on the evolutions of a single static black 
hole, of a single moving black hole with and without spin, 
and of the head-on collision of a BBH, with and without 
a massless scalar field. We summarize and discuss the 
implications of our findings in Sec. IVII 



II. FORMULATION 

A. Review of the Basic Equations 

The code is based on the BSSN formalism 10], which 
is a conformal-traceless "3-1-1" formulation of the Ein- 
stein equations. In this formalism, the spacetime is 
decomposed into three-dimensional spacelike slices, de- 
scribed by a three-metric 7^ ; its embedding in the four- 
dimensional spacetime is specified by the extrinsic curva- 
ture Kij and the variables, the lapse a and shift vector 
/3', that specify a coordinate system. Our conventions 
are that Latin indices run over 1, 2, 3, whereas Greek in- 
dices run over 0, 1, 2, 3. Throughout the paper we adopt 
geometrical units with G = c = 1. In this paper we fol- 
low the notations of [2^. The metric 7.^ is conformally 
transformed via 

= -ji ln7, 7ij = e^*'''jij, (1) 



where 7 denotes the determinant of the metric jij . The 
conformal exponent (jj is evolved as an independent vari- 
able, whereas 'jij is subjected to the constraint that the 
determinant of 7^ is chosen to be unimodular, i.e., 7=1. 
The extrinsic curvature is subjected to the same confor- 
mal transformation, and its trace K is evolved as an in- 
dependent variable. That is, in place of Kij we evolve: 

K = j'^K,j, A,,=e-^'^K,,-^^,,K. (2) 

Similar to the conformal metric, Aij is subjected to a 
constraint that Aij is traceless, i.e., tiAij = 0. New 
evolution variables, i.e., the conformal connections, 

r^-f\,, (3) 

are introduced, defined in terms of the contraction of the 
spatial derivative of the inverse conformal three-metric 



With these dynamical variables the evolution equa- 
tions read 
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Here p, s, Si, Sij are source terms which come from mat- 
ter. For a vacuum spacetime p = s = Si = Sij =0. In the 
above evolution equations Di is the covariant derivative 
associated with the three-metric 7^ , and "TF" indicates 
the trace-free part of tensor objects. The Ricci tensor 
Rij is given as 

i?,, = (9) 

+7™"(2r''m(irj)fc„ + T''inTkmj), (10) 

+4D,(j,D,^-Aj,,D^cbDk(f>. (11) 

The Einstein equations also lead to a set of physical con- 
straint equations that are satisfied within each spacelike 
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slice: 

e-'^'^{R - 8&D^<I> - 8D''<f>Dk(t)) 

+ ^K^ - AijA'^ - Wirp = 0, (12) 
_ > > . 2 ~ 

D'Aij + 6AijD'(t) - -DjK - Stts^ = 0, (13) 
o 

which are usually referred to as the Hamiltonian and the 
momentum constraints. Here R = R^i is the conformal 
Ricci scalar on a three-dimensional time slice, and Di 
is the covariant derivative associated with the conformal 
three-metric . Besides being used to obtain the evolu- 
tion equations ([6]) and ([8]) in the BSSN formulation, the 
Hamiltonian and the momentum constraints are also ap- 
plied to the volume integrals of the ADM mass and the 
angular momentum, respectively [26| : 

^ = / - 8r^a,e"*)dS, (14) 

167r /go 

= / A^x[e'^{lQ^p + A,,A^= -\k^) 

+f^fc-e*i?], (15) 
J^ = ^e./ i e^-^x^i^dE^ (16) 

-^x^A£™f™,fc + 87ra;^'sfe)], (17) 

where dEi = (l/2)eijfcdxMx'^. These two global quanti- 
ties are useful tools for the system diagnostics to validate 
the calculations. The volume integral (fT5|) is slightly dif- 
ferent from the one in ^26] due to the further application 
of the unimodular determinant of the conformal metric 
(fT8)) . Refer to Appendix [Cl for the details. 

B. Equation Adjustments 

The specific choice of evolution variables introduces 
five additional constraints, 

7-1 = 0, (18) 
tri.j = 0, (19) 
f'+f\j = 0. (20) 

Our code actively enforces the algebraic constraints 
and (|19p by replacing "fij and Aij with the following: 

% ^ (21) 

Aij > Aij — —jijtrAij. (22) 

To enforce Eq. all the undifferentiated in the evo- 
lution equations are substituted with —7*^ j . 



As to the variable choice for the conformal factor, 
the alternative x, first proposed in Q, has been widely 
adopted. In the x method the conformal exponent 4> 
[which is O(lnr) near the puncture] is replaced with a 
new variable x = e^^"^ [which is 0(r*) near the punc- 
ture] . X grows linearly near the puncture during the time 
evolution; such linear behavior leads to a more accurate 
evolution near the puncture. In the x method, equation 
dH) is replaced by 

dtX^lx{aK-(3\,)+(3\,,. (23) 

Note that i^,, = -x,i/4x and (j)^ij = X,jXj74x^ - X.ij/'^X 
are applied to the evolution equations ([6]), ^ [via 
Eq. pT|) ]. and ([5]). In these substitutions, the divisions 
by X need to be taken care of in the numerical implemen- 
tation to avoid division by zero or unphysically negative 
values of X- In HI] a small e is set to replace x in division 
if X is less than e. In [2II W = e~^'^ is chosen to be the 
conformal factor variable instead of x, to avoid the effect 
of unphysical negative values of x on the evolution of the 
other variables [53]. However, we did not encounter any 
such difhculty in the work for this paper, therefore it is 
not necessary to apply the aforementioned modifications, 
although we anticipate the appearance of this difficulty 
in some complicated scenarios in future work. 

C. Gauge Conditions 

As mentioned in the introduction, the gauge conditions 
are important for the numerical simulations of dynami- 
cal spacetime, and this is especially true for the moving 
puncture method. The Bona-Masso type slicing gauge 
conditions [21I ] for the lapse function and many driver 
gauge c onditions (e.g., the F-driver) for the shift vector 
|2a . [23] are currently the main type of gauge conditions 
used in the punctured black hole calculations. In this 
work, wc will only focus on these types of gauge condi- 
tions, which can be written as 

dta = ~2aK + XlP'a.^, (24) 

dtP' = ^f{a)B' + \2(3'(3\,, (25) 

dtB' = dtr-r^B' +X3f3^B\,-X^p^f\,, (26) 

where rj and the four A's are the parameters to be chosen, 
and f{a) is a function of a. Here the A's can only be set 
to be or 1. We set f{a) = 1 in all the cases except in 
Sec. IV Al where f{a) = a. The gauge conditions used 
for moving black holes in the literature include: (1) Ai ~ 
1, A2 A3 A4 = (e^ see [ig^); and (2) Ai = A2 = 
A3 = A4 = 1 (e.g., see [ill, [l3|), with the proper ry's. In 
[22| . the authors investigated several cases of the above 
gauge equations. In [33|, the authors discussed the above 
case (2) analytically. In this work, we will explore this 
problem more thoroughly (see the following sections for 
details). In particular, we are concerned about the effect 



of the advection terms on the stabihty and accuracy of 
the evolution, and we try to find the viable gauges for 
moving and/or spinning black holes. Throughout this 
paper we will fix 77 = 2. 



III. NUMERICAL METHOD 

In this section, we briefly describe the key numeri- 
cal techniques used in this work. For the discretization, 
our code uses the cell-centered method, which takes the 
data to be defined at the center of the spatial grid cell. 
We also use a finite centered-differencing method with 
fourth-order accuracy to approximate the spatial deriva- 
tives, which closely follows [2^. In the temporal part, we 
use the iterative Crank-Nicholson method for time inte- 
gration, which gives a second-order accuracy [s^l- We 
take di=(Courant factor) xdx, and the Courant factor is 
set to be 1/4. 

We apply the standard centered finite differencing ap- 
proximation to all spatial derivatives except the advec- 
tion terms (i.e., the terms of the form f3^djF). For these 
advection terms we use the following fourth-order lop- 
sided stencils: 

+10F,j^fe-|-3F,+ij-fc) for /3"<0, (27) 

dxFi,j,k — Y2da; ~ ^ 18-Fi+ij,fc 

-10^,,,- fc - 3F,_i,,, fc) for r>0, (28) 

along the s-direction. The stencils are similar along the 
y- and z-dircctions. We also install in the code a Kreiss- 
Oliger dissipation [31] of the form 

dtF -> RHS + e(-l)"/2 h'l+'D'^l^+^Df^+^F, (29) 

where RHS represents the corresponding evolution equa- 
tion for F, hi is the grid spacing in the ith direction, 
Di^ and are the forward and backward differenc- 
ing operators in the ith-direction, n is the order of the 
finite difference used to evaluate the RHS, and e is the 
dissipation coefficient to be chosen in various cases. 

In order to increase the numerical resolution without 
increasing the computation cost much, the mesh refine- 
ment method is used in the numerical simulations. We 
use the GrACE package [IJl to implement both the mesh 
refinement and the parallelization in our code. This pack- 
age is able to deal with the adaptive system, consider- 
ing both partitioning and load-balancing for distributed 
adaptive mesh refinement applications. However, we only 
use fixed mesh refinement (FMR) in this work. The com- 
putational domain is represented by a hierarchy of nested 
Cartesian grids; we adopt the cell-centered scheme, and 
the grid hierarchy follows the Berger and Oliger algo- 
rithm [131 . The hierarchy consists of L levels of refine- 
ment indexed by = 0, ... , L — 1, for which ^ = is 



4 



— «K 


— «K 


— IK 




* 




— «K 










— «K 




♦ 
♦ 






— <K 



















FIG. 1: (closely follows Fig.l of [S^) Two-dimensional dia- 
grams for the guard cell filling. The thick vertical lines in 
both panels represent the refinement boundaries separating 
fine and coarse grid regions. The left panel shows the first 
step, in which one of the coarse grid cells (red circles with 
black rim) are filled using a quartic interpolation across 25 in- 
terior fine grid cells (green diamonds). The right panel shows 
the second step in which two fine grid guard cells (red dia- 
monds with black rim) are filled using quartic interpolations 
across 25 coarse grid values (circles). These coarse grid val- 
ues include two layers of guard cells (green circles), obtained 
from the coarse grid region to the right of the interface, and 
three layers of interior cells (green circles with black rims). 
The final step (not shown in the figure) is to use "derivative 
matching" to fill the guard cells for the coarse grid. 



the coarsest level and £ = L — 1 is the finest one. A 
refinement level consists of one or more Cartesian grids 
with constant grid spacing he on level £. All grid blocks 
have the same logical structure, and refinements are bi- 
sected in each coordinate direction. A refinement factor 
of 2 is used such that h£ — ho/ 2^. The grids are properly 
nested such that the coordinate extent of any grid at level 
£, £ > 0, is completely covered by the grids at level £—1. 
We do not refine the time step. So the time step for each 
level is set as dt = (Courant factor) x /imin, where /imin 
is the resolution of the finest level. Therefore, there is no 
interpolation of data between different time slices. 

Another problem in the process of the mesh refine- 
ment is how to treat the refinement boundary well to 
avoid possible numerical noise. Here we follow the guard 
cell scheme described in |33|. This method has three 
steps, shown in Fig. [TJ In the first step, interior fine 
grid cells are used to fill the interior grid cells of the 
next lower level. This restriction operation is depicted 
for the case of two spatial dimensions in the left panel 
of Fig.[TJ The restriction is basically a three-dimensional 
interpolation in this work, and is accurate to fourth-order 
in the grid spacing to match the accuracy of the finite- 
differencing scheme. As in the left panel of Fig. [Tl the 
values of the coarser grid cells (red circle with black rim) 
located within the finer grid cells (green diamonds) are 
filled with the finer values by quartic interpolations. The 
stencil of the finer cells needs to take care to ensure that 
only interior fine grid points, and no fine grid guard cells, 
are used in this first step. Secondly, the fine grid guard 
cells (red diamonds with black rims) are filled by pro- 
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longation from the grid of the next lower level. Before 
the prolongation, the coarser grid updates its own cells 
(green circles with black rim in the right panel of Fig. [1]) 
from the finer grids in the first step. The stencil used 
in the prolongation operation is shown in the right panel 
of Fig. [TJ The prolongation is a three-dimensional in- 
terpolation, and is also accurate to fourth-order, like the 
restriction in the first step. In this case, the coarser grid 
stencil includes two layers of guard cells (green circles), 
as well as its updated interior grid points (green circles 
with black rims). In the last step, the coarser guard cells 
close to the interface (bold line) are filled by using deriva- 
tive matching, the difference between the finer cell and 
the neighbor finer guard cell across the interface matches 
the difference between the coarser cell and the neighbor 
coarser cell across the interface 1541. 



with 



U2 



U3 



4062 



(15 + 117£-79£2-|-43^3 



20 



-14£* + 2£5 + 841n^/6), (33) 



U4 = —{l0-25i + 2lf-6£^), 

where b — 2r/m, £ — 1/(1-1- b), jip = P ■ n/P and 
/i5 — S ■ n/S. For the approximate solution (|32p of a 
moving black hole with spin, the ADM mass, the linear 
momentum, and the angular momentum are 



IV. INITIAL DATA FOR PUNCTURED BLACK 
HOLES 



For the initial data of punctured black holes, we con- 
sider the Bowen-York type initial data, in which the max- 
imal slicing and the conformal flat form are adopted [l^ . 
Let ip be the conformal factor, = e"^. The conformal 
extrinsic curvature reads 



A, 



V.J - nln])Pfni + -^[,6,)^, 5^^^1,(30) 



5 P2 2 52 

Madm = "1+77 ^ 

8 m 5 m'^ 

Pabm = P, 
Sadm — S. 



(34) 

(35) 
(36) 



On the other hand, when P = while S is very large, 
the conformal factor can be approximated as (see Ap- 
pendix]^ for more detail) [s^ 



(6^^) 



2U/8 



(37) 



In this paper we will fix m = 1 for all single black hole 
simulations. 



where fij is the three flat metric; P} and Sj are constant 
vectors, standing for the linear momentum and spin mo- 
mentum of the /-th black hole respectively; Uj is the ra- 
dial normal vector with respect to fij , which points from 
the position of the /-th black hole to the space point. In 
the puncture method described in [33| ,'0 = 1 + S^+" 
with mass parameter mj for the /-th black hole, and u 
is determined by 

{dl + dl + dl)u = -\k''k,A^ + E ^ + (31) 

If all black holes are at rest and spinless, ([30]) implies 
Kij = 0, so u = 0. In the head-on collision case, we 
will implement this kind of initial data with mi = 1712 = 
0.5. For a single black hole, when the linear momentum 
and the spin are small, we can solve the above equation 
approximately as (see Appendix El for more detail) [s^ 

U = —[Ui ^ U2{i^il - I)] + 6^^2(1 +//|) 

TO* 

+ ^PxS-n, (32) 



V. NUMERICAL RESULTS 

In this section we report the numerical results for: (1) 
a single moving black hole without and with spin: The 
moving action and the spinning action of a single black 
hole are fundamental elements for BBH simulations. As 
our code aims at simulating BBH coalescence, evolving 
a single moving and spinning black hole becomes an es- 
sential test. In addition, since the gauge choice is crit- 
ical for the moving puncture method, we would like to 
study if there are any other gauge conditions which can 
also support the moving puncture technique, besides the 
known ones. Our main achievement in this part of the 
work is that we discover one new gauge condition be- 
sides the known ones which can support moving punc- 
ture black hole simulations. The results on the gauge 
condition tests are listed in Table [H where Gauge VII 
is the aforementioned new set of gauge conditions. The 
successes of the gauge usage in other groups' work, e.g., 
[ilSIiililEii, are also reconfirmed in Table |T1 (2) 
the head-on collisions of BBHs: The head-on collision of 
a BBH system is the simplest dynamical spacetime in 
which a complete gravitational waveform of the merger 
of two black holes could be produced. Therefore, we use 
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TABLE I: The gauge choices tested in this work correspond to 
Eqs. (gH), (123 and ((251); "V" stands for the gauge condition 
with the corresponding advection term, while "x" stands for 
without. 
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FIG. 2: The root mean square of the change in the trace 
of extrinsic curvature between consecutive time steps as a 
function of time in the static case with equatorial symmetry. 
The solid (red) line is the result with the conventional setting. 
The dashed line is the result with the modifications suggested 
in [2^. This shows that the both settings give stable and 
convergent results. 



this scenario to examine the performance of the code. 
Besides, in order to go beyond the cases of a head-on col- 
lision in vacuum, the case of a head-on collision perturbed 
by a massless scalar field is also studied. With such kinds 
of cases we try to understand qualitatively the effect of 
the existence of neutral matter on the collision of a BBH, 
especially on its gravitational radiation. It is shown that 
the waveform could be affected significantly by the scalar 
field. 



A. Static Black Hole 

Although the BSSN formulation with the "1-hlog" 
lapse condition and the F-driver shift condition has been 
shown to be well-posed and hyperbolic [13, , it is still 
useful to confirm the stability and convergence of the for- 
mulation and the applied modifications before we move 
forwards to the moving/spinning black hole cases in the 
following subsections. Therefore, for the equation ad- 
justments, we enforce the constraints p ^ - (|^D)) by us- 
ing Eqs. (PTjl and (P^ . and the substitution of the con- 
formal connection with — 7*-' j, as described in Sec. 
IIIBI For the gauge condition, Eqs. P^ - (PS)) are ap- 
plied with /(a) = a and the parameter choice A3 = 0, 
Ai = A2 = A4 = 1, which is close to Gauge VII in Ta- 
ble |T1 i.e., the newly viable gauge condition (see Section 
IV Bp . The grid width h — 0.2 and the outer boundary 
is ±16, ±16, 16, respectively, assuming equatorial sym- 
metry. In this simple case we only consider a unigrid for 
the computational domain. 

The result for such a conventional setting is shown in 



Fig. O This figure shows a log plot for the root mean 
square (r.m.s.) of the changes in the trace of extrinsic 
curvature K (the solid red line) between consecutive time 
steps. In the plot the curve of the change in K rises 
around t = 200 during the period of settlement. The 
change in K decreases exponentially afterwards, without 
a sign of rise to the end of the run. This indicates that the 
conventional settings give a stable evolution for a single 
static black hole. The result is also consistent with the 
analytic understanding of the BSSN formulation and of 
the gauge choice. 

Meanwhile, some modifications [2§|, especially the en- 
forcement of the constraints (|18P " (pn|) . have been shown 
to be at least as good as the conventional ones in the 
numerical result [l3|. We are therefore interested in un- 
derstanding the performance of the code with the modi- 
fications, and the comparison between these two sets. 

Briefly, the modifications in [26l | are summarized as fol- 
lows: Instead of treating all components of 7ij equally, 
only five of the six components of are evolved dynam- 
ically, and the zz-component is computed using Eq. ([18]), 



= 1 



Ivvllz 



Ixylyzlxz + Ixxl"^ 



Ixxlyy 



Ixy 



(38) 



Similarly, A^z is determined from the other five compo- 
nents of Aij using Eq. (fTT 



Azz 



;^zz 



(39) 



Instead of substituting for the undifferentiated conformal 
connection F' by — T"*-* j according to the constraint (|20p . 
the constraint is added to the evolution equation ([5]) of 
via 



d.r = rhs of 



(e + 2/3)(r+f^,)/5',fe, (40) 



where ^ is usually chosen to be 2/3. 

With otherwise the same settings as in the conven- 
tional ones, the result with the modifications in [26j is 
also indicated by the dashed line in Fig. [21 Without too 



7 



long a settlement period, the change in K decreases expo- 
nentially to the end of the run. A comparison of these two 
sets of modifications shows that they are equally good in 
converging the evolution to a numerically stable state, 
although the modifications in [HI give a better settle- 
ment at the early stage of the evolution. This indicates 
that there is still room for modifying the BSSN formu- 
lation to achieve better stability and convergence. The 
major purpose of this work is to build a reliable code for 
BBH simulations, therefore we will stick to the conven- 
tional modifications in the rest of this paper, although 
the modifications in [2^ might show some subtle advan- 
tage in stabilization over the conventional ones. 

B. Moving Black Hole without Spin 

We now study the cases of a moving black hole without 
spin, which are similar to the ones in (22.] . Three levels 
of grids are used in this and the next two subsections. 
The outer boundary of the coarsest level is set at ±16, 
and the boundaries of the finer levels are located at ±8 
and ±4, respectively. The gridwidth for the highest level 
is 1/8. The black hole is located at (0,0,-3) initially 
with the linear momentum vector P = (0,0,1). The 
gauge choices tested in this work are listed in Table H] 
Gauges I, II, IV, V, VI and IX (Figs. 3, 5, 7, 8, 6, 10 in 
, respectively) have been tested in |22il . Differing from 
the numerical initial data used in [22], the approximate 
analytic initial data described in Sec. IIVI is used in the 
tests. Nevertheless, our results are consistent with the 
results in pi ]. 

We summarize our results obtained for the nine tested 
gauges in Fig. O In each panel of Fig. [31 the number 
on the upper-right corner indicates the case with the 
same gauge number in Table [H and the result of that 
case is plotted in the panel. In each case, the black hole 
moves along the z-axis, and the conformal connection 
^-component T^ (the red solid line), the lapse function 
a (the dashed line), the conformal exponent (the ma- 
genta dot-dashed line) , and the shift vector z-component 
f3^ (the blue dot-dot-dashed line) are chosen as the mon- 
itors for the stability of each run. The profiles of these 
variables are recorded at time i = 30 in each panel, ex- 
cept for Panel VI, where the recorded time is t — 12. In 
Panel I, it can be seen that a, (p, and f3^ all behave well, 
but has some tail of ripple behind the black hole. The 
ripple implies a rising instability in the evolution. It is 
known that adding the a advection term in Eq. (|24p . i.e., 
Ai = 1, suppresses this unstable mode. We verify it by 
comparing the profiles of the variables in Panel I and II. 
In Panel II, all variables behave well and there is no rip- 
ple of noise for the curve of F^. This result is consistent 
with the one in [si]. We then set Ai = 1 in Eq. ((24|) in 
the following cases, since it is necessary for the stability 
of a moving black hole. 

In Panel III, it is shown that the addition of the advec- 
tion term of /? in Eq. ([^5)1 . i.e., Ai = A2 = 1, will result 



in a big bump behind the black hole. This instability is 
somehow strong, such that it causes the ill-behavior of a 
and (f). Adding the advection term of B in Eq. P^ . i.e., 
Ai = A3 = 1, will bring in a small tail, as seen in Panel 
IV, although such an instability does not seemingly af- 
fect a, (j), and /3^, for they behave well in the plot. The 
subtraction of the advection term of F' in Eq. ([^5)1 . i.e., 
Ai = A4 = 1, will result in a "distorted" profile for F^ 
in Panel V. (A gauge choice close to Gauge V has been 
used to simulate the inspiral of a BBH system with a 
small initial separation in 9].) From the above three 
cases, we see that solely adding the advection term of /3 
or B, or subtracting the advection term of F* will bring 
in instability in general. This understanding leads us to 
try combinations of these three cases. 

Panel VI shows that the combination of the /3 and the 
B advection term additions, i.e., Ai = A2 = A3 = 1, re- 
sults in a profile with high frequency noise for F^ around 
the black hole and makes an even stronger instability. 
The code crashes soon after time t = 12 when the black 
hole has only moved slightly. However, in Panel VII, 
it is shown that the combination of the /3 advection 
term addition and F* advection term subtraction, i.e., 
Ai = A2 = A4 = 1, can suppress the unstable modes in- 
troduced by the sole usage of these two terms (see Gauges 
HI and V) . We can see that the curve profiles of all the 
variables in Panel VII look almost the same as those in 
Panel II. It is interesting to note that, according to our 
literature search. Gauge VII has never been previously 
used in any BBH simulations. Therefore, it deserves fur- 
ther study. In Panel VIII, we consider the combination 
of adding the B advection term and subtracting the F' 
advection term, i.e., Ai = A3 = A4 = 1. The performance 
in Panel VIII shows a set of curve profiles which are sim- 
ilar to those obtained with Gauge V. Finally, we consider 
in Panel IX the combination of all three advection terms 
in the shift equation, i.e., Ai = A2 = A3 = A4 = 1. 
The combination of the BSSN equations with Gauge IX 
has been proven to be strongly hyperbolic in the sense 
of first-order in time, second-order in space systems j37| . 
and thus yields a well-posed initial-value problem. Panel 
IX shows that all the variables behave very well and 
smoothly. The curve profiles in Panels II, VII, and IX 
look almost the same. 

The linearized analysis of the BSSN formulation with 
the gauge conditions in [23] shows that both Gauges II 
and VII have zero-speed modes. However, from the re- 
sults described in this section, we can not distinguish the 
difference between Gauges II, VII and IX. Furthermore, 
Fig. 5 of [l^, which corresponds to gauge II, shows no 
zero-speed modes at all. We conjecture that the nonlin- 
earity of the full theory could eliminate the zero-speed 
modes for Gauges II and VII in general. In fact, Gauge 
II [H, as well as Gauge IX [H [H, [H, has been 
successfully used in the simulations of black hole evolu- 
tion. And Gauge VII is as good as Gauges II and IX 
for black hole simulations, at least in the cases tested in 
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FIG. 3: The gauge tests for a moving black hole without spin (velocity v « 0.615). The profiles of several dynamical variables 
at t = 30 are shown, except for case VI, for which the time is t = 12. The horizontal axis is the z-axis, the moving direction 
of the black hole. The vertical axis is the corresponding value for different variables: the solid (red) line is F^; the dashed 
line is a — 1; the (magenta) dot-dashed line is <j); the (blue) dot-dot-dashed line is (i^ . The different panels correspond to the 
different cases in which the gauge number is marked on the upper-right corner of the panel and they are also listed in Table U 
For Gauges III, V, and VIII, the results show some ill behavior explicitly while Gauges I and IV have tails of noise behind the 
black hole. The rest of the three gauge choices II, VII, and IX give almost the same well-behaved result. The results of Gauges 
II and IX are consistent with other research groups' work. Gauge VII is found to work well with the moving puncture method 
in this work. 
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FIG. 4: The gauge tests for a moving black hole with spin parallel (upper panel) and perpendicular (lower panel) to the moving 
direction (velocity v ~ 0.5, the specific angular momentum a « 0.5f02). The profiles of several dynamical variables at t = 30 
are shown as in Fig. [3] From left to right, the figures correspond to the results with Gauges II, VII and IX listed in Table U 
respectively. These three "good" gauges give almost the same results. 



this work. Therefore, we can also expect that Gauge VII 
is very likely to be viable in the generic cases of BBH 
evolution. 

Someone may raise the doubt that the moving velocity 
is not large enough to excite the zero speed mode for the 
newly found gauge condition here. This is not the case 
in fact, at least for the approximate initial data. Consid- 
ering the ADM mass (j34[) . we have the maximal moving 
velocity v = Pabm/Mabm for a black hole without spin 
when P — 2\/lO/5. And the moving velocity {v « 0.615) 
for the above tested moving black hole with the linear 
momentum P = I almost equals this maximal velocity 
(v w 0.63). We have tested this maximal moving veloc- 
ity also. The result gives the same conclusion mentioned 
above. Someone may also raise the doubt that the zero 
speed modes for Gauge VII mentioned in [22] might not 
be excited in this tested scenario. Therefore, we test it 
further, with both the case of a moving black hole with 
spin and the case of a high-spin black hole, in the follow- 
ing subsections. Our results will show that Gauge VII, 
as well as Gauges II and IX, passes these two tests. 



C. Moving Black Hole with Spin 

In this subsection, we investigate the effect of spin on 
a moving black hole. Firstly, we set the spin direction 
of the black hole to be parallel to its moving direction, 
specifically, P = (0,0,1) and S = (0,0,1). Therefore, 
the only difference between the setting here and in the 
previous subsection is that, in this subsection, the black 
hole has a spin of amplitude 1 which is along the moving 
direction. Here we only test the "good" gauge choices, 
i.e.. Gauges II, VII and IX. The curve profiles at i = 30 
are presented in the upper panels of Fig. 3) We find 
that, with any of these three gauges, the code can sta- 
bly simulate a spinning black hole. The curves in the 
three upper panels look the same, not showing a sign 
of instability. Compared with Fig. [31 we find that the 
black holes in this case move slower {v « 0.5) than the 
spinless ones {v « 0.615) with the same gauge choices 
described in the previous subsection. This is consistent 
with the theoretical prediction: From Eq. (|34p we can 
see that the moving black hole with a nonvanishing spin 
will have a larger ADM mass than the spinless one due 
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FIG. 5: The gauge tests for a single rapid- rotating black hole (spin parameter a ~ 0.9). The comparison of profiles of the three 
"good" gauges at t = 30 is shown. From top to bottom, and from left to right, the panels show the conformal factor (j), the 
lapse a — 1, the conformal connection F^ and the shift respectively. These three gauge choices give almost the same results. 



to the rotational energy. For the moving velocity being 
V = Padm/Mabm, a spinning black hole will move slower 
than a spinless one with the same linear momentum. Our 
conclusion in this test is that Gauges II, VII and IX all 
handle this scenario well, since these three gauges give 
almost the same result. 

Secondly, we set the spin direction of the black hole 
to be perpendicular to the moving direction, specifically, 
P = (0,0,1) and S = (1,0,0). Similarly, we only test 
Gauges II, VII and IX. The curve profiles at i = 30 are 
presented in the lower panels of Fig. [D The results are 
similar to those obtained in the parallel-spin case. From 
Eq. the moving velocity of the black hole for this 

case is the same as the velocity in the parallel-spin case, 
i.e., V « 0.5. We can read this from the results presented 
in Fig. m In this case, the peak amplitude of the variable 
(j) becomes smaller than in the parallel-spin case. How- 
ever, this difference basically comes from the spin ori- 
entation; it does not bring any instability to the results 
with the three "good" gauges II, VII and IX. These three 
"good" gauges give essentially the same stable behavior. 

Similar to the moving velocity being bounded, the 
magnitude of the specific angular momentum, a = 



S/M'pjyj^, has a maximal value for the approximate ini- 
tial data ([5^ . For S — \/5/6, the black hole has a 
maximal magnitude a — 3-\/30/32 « 0.5135. The mag- 
nitude of a in both the parallel-spin and perpendicular- 
spin cases is a « 0.5102. Therefore, the magnitude of a 
in these two cases is very close to the maximal magni- 
tude of the specific angular momentum that the initial 
data can have. In the next subsection we would like to 
test these three gauge choices with an even higher spin 
magnitude of a. 



D. Rapidly-rotating black hole 

We have tested Gauge VII in the spinless and the 
spinning moving black hole cases. We find that, with 
this gauge choice, the code can simulate these dynam- 
ical spacetimes well. Since experience tells us that a 
higher spin is more likely to trigger a quicker instability 
in the numerical simulations, one might doubt whether 
the newly found gauge condition can handle a high-spin 
case well. In this subsection, we test the three "good" 
gauge choices with a rapidly-rotating black hole. The ini- 
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tial data used is another approximate analytic solution 
for the Bowen-York initial data. We describe the details 
in Appendix [Al This initial data is similar to that used 
in Sec. IV C[ except that the conformal factor tp is given 
by Eq. ^ rather than Eq. ([321) (sl, and F = 0. 

Here we set the angular momentum vector to be 5 = 
(0,0,10000), which results in the specific angular mo- 
mentum a « 0.9 |36|]. In this case, we again only care 
about the corresponding behavior of the runs with the 
three "good" gauge choices. The outer boundary is set 
at r = 128, and six levels of grids for FMR are used in 
the runs. The curve profiles at t = 30 are presented in 
Fig. m These three "good" gauge choices can all give a 
stable and accurate simulation for this highly spinning 
single black hole case. Looking at the plots in this figure, 
it is difficult to distinguish between the results obtained 
with these three gauges. The curves overlap each other 
well. The profiles of and are consistent with the 
experience that they almost give the same shape. Due 
to the limits on computational resources, the runs are 
stopped at t = 30. This might raise the doubt whether 
the runtime is long enough to excite an unstable mode. 
Since the spin of the black hole in this case is close to the 
maximal one (cmax ~ 0.9282) that a punctured black hole 
can have in a Bowen-York type initial data [36] , we expect 
that some difference between these runs, caused from in- 
stability, will appear in this period of time. However, we 
find that the results with these three gauge choices are 
almost the same. This should indicate that Gauge VH is 
as good as the other two in this case, unless all three of 
these gauge choices cause the same instability, which is 
highly unlikely. Meanwhile, the results in Sec. lV Al can be 
regarded as complementary to the one in this subsection 
(about the the long-term stability issue of Gauge VII). 
Therefore we conclude that the new gauge choice, Gauge 
VII, survives in this high-spin test. 



E. Head-on Collision of two equal mass black holes 

In the previous subsections, we have presented the nu- 
merical simulations for a single black hole with/without 
spin with the different gauge choices. In order to test our 
code as well as the three gauge choices discussed above 
further, we present in this and the next subsections the 
numerical results about the head-on collisions of a BBH 
system, which is the simplest case for BBHs. We use 
time-symmetric initial data for the two black holes, i.e., 
the so-called Brill-Lindquist initial data [4l| . Specifically, 
the initial data takes the form 

2|r-ci| 2|r-C2| 
7u = K^O, =0. (42) 

where mi and m2 are the mass parameters for the two 
black holes, ci and C2 are the positions of the black holes, 
and fij stands for the flat three- metric. In this case. 
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FIG. 6: The (^=2, m=0) mode of r^4 extracted from a head- 
on collision of the Brill-Lindquist initial data starting from 
the approximate ISCO separation d = 2.303 at four different 
radii r=20, 30, 40 and 50. Considering the velocity of the 
gravitational wave, the time delay has been shifted in the 
plot. 



we set the mass parameters to be mi = m2 = 0.5 and 
the positions of the two black holes at (0, 0, ±1.1515), 
respectively. This value has been used in [ij, IH, and 
corresponds to an initial separation of the black holes 
equal to that of an approximate ISCO configuration. In 
other words, our initial data corresponds to two identical 
black holes which have no spin and no linear momentum, 
attracting each other from rest at the ISCO. 

In this and the next subsections, the computational 
domain is ±64 x ±64 x 64 and 64 x 64 x 32 grid points 
are used in every level. To reduce the computational load 
equatorial symmetry is assumed. Six levels of grids are 
used for the mesh refinement. The refinement boundaries 
are placed at 32, 16, 8, 4, and 2. As mentioned in Sec lIIBl 
the version of the evolution equations are also available 
in our code. In this and the next subsections, we test 
both the 0- version and the X" version of the code in the 
head-on collision cases. The results obtained from these 
two versions are consistent in the scenarios fssl] . 

Although conceptually simple, the head-on collision 
of a BBH system is still a highly dynamical process in 
spacetime. During the collision, the system will emit a 
complete gravitational waveform from the merger of two 
black holes. To quantify the gravitational waveform, we 
use the Newman- Penrose scalar ^"4. The method of com- 
puting this quantity is described in Appendix [Bl Due to 
the symmetric property of this BBH system, the gravi- 
tational wave has only the (i?=2,m=0) mode. First, we 
consider Gauge IX, which has been widely adopted in 
simulations of BBH systems. We extract the waveform 
at r=20, 30, 40, and 50. Since the leading order of ^'4 
is 1/r asymptotically, r^'4 should be independent of the 
extraction position, except for the time delay due to the 
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FIG. 7: The differences between the (£=2,m=0) mode of the 
waveform for different gauge choices extracted at r = 30. 
The sohd hne is the difference between Gauges II and IX, the 
(red) dashed line is the difference between Gauges VII and IX, 
and the (ohve) dot-dot-dashed hne is the difference between 
Gauges II and VII. The largest difference is roughly 0.2% of 
the amplitude of 



FIG. 8: The (^=2,m=0) mode of the waveform r»I'4 for a 
massless scalar field perturbing the head-on collision of two 
identical black holes without spin or linear momentum. A 
is the strength of the perturbation. A — 0.0 means without 
perturbation. The extraction radius is r = 30. The time 
delay effect of the perturbation is clear. The perturbation 
will efficiently amplify the waveform. 



gravitational wave propagation. We consider the velocity 
of the gravitational wave to be the speed of light and thus 
subtract the time delay accordingly. The waveforms are 
plotted in Fig. [S] The result is quantitatively consistent 
[56.] with the result reported in [l^, . Initially there 
seems to be some small amplitude oscillations before the 
larger oscillations for the monitors at larger radii r, but 
not for the one at smaller r. This noise mainly comes 
from the reduced accuracy of the evolution in the coarser 
grids. 

We next study the effect of the gauge conditions dis- 
cussed in the previous subsections in the head-on colli- 
sion. The differences of the waveforms with the gauge 
choices are shown in Fig. [T] As expected, the codes with 
these three "good" gauge conditions can all evolve the 
head-on collision process stably. However, the gauge con- 
ditions indeed affect the waveform [45| . The highest peak 
of the difference between Gauges II and IX, A(r\l/4)ii_ix, 
is about 4 x 10"^. From the plot, we can see that the 
pattern and the amplitude of A(r5'4)vii-ix, the (red) 
dashed line, is close to A(r4'4)ii_ix, the solid line. The 
amplitude of the wave is about 0.02, so the relative dif- 
ference is about 0.2% of the amplitude of But it is 
interesting to note that the amplitude of A(r\E'4)ii_vii, 
the (green) dot-dashed line, is smaller, only about 0.05%. 
We note that these differences are smaller than the dif- 
ference of wave form resulting from the different initial 
lapse profile (for example a = 1, l/tp'^, etc.) which 

is typically 1% of the amplitude of rvE'4. 



F. Head-on Collision Perturbed by a Massless 
Scalar Field 

In the previous case, we studied the head-on collision 
of a BBH system without matter. However, matter is at- 
tracted into strong gravitational systems. Thus matter 
also plays an important role in binary compact objects. 
The evolution of a dynamic spacetime generally needs 
to include matter. On the other hand, most of the as- 
trophysical systems, including BBH mergers, are usually 
surrounded by an accretion disk or some kind of matter. 
Therefore, it is interesting to ask how this matter will 
affect the gravitational wave signal, which is expected to 
be detected in the near future. In order to check the 
consistency with the matter coupling in the code and to 
investigate the effect on the gravitational waveform by 
matter, we study the head-on collision "perturbed" by a 
scalar field in this subsection. Here "perturbed" means 
that the amplitude of the scalar field is small, and we 
do not need to solve the constraint equation to get exact 
initial data. Instead, we set the dynamical variables of 
geometry to be identical to those in the previous sub- 
section. Furthermore, the matter part of the dynamical 
variables are set independently. However when we evolve 
this initial data, we do not do any approximation, that is 
to say we solve the fully coupled Einstein-Klein-Golden 
equation numerically. 

The evolution equation of a scalar field is described in 
Appendix [D1 For each simulation, we provide the follow- 
ing initial data for $ and 9t$: the profiles of are a rest 
sphere shaped scalar field located at the center of the two 
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black holes 

$(t = 0) = Ae-''', 9t$(t = 0)==0, (43) 

where A is the amplitude of the scalar field. We test 
the scalar field with different amplitudes to see its effect 
on the waveform. The amplitudes are set to be A=0.02, 
0.04, 0.06, and 0.08. The results are plotted in Fig. El 
It is clear to see that the waveform is delayed with the 
existence of the scalar field, and the larger the amplitude 
of the scalar field is, the longer the delay time. In the 
meantime, the amplitude of the waveform also becomes 
larger. This phenomenon can be understood as follows: 
The scalar field is located at the middle of the black hole. 
Some part of the scalar field will escape outwards when 
it evolves. This escaping part of the scalar field will de- 
lay the motion of the two black holes toward each other. 
Meanwhile, some part of the scalar field will be absorbed 
into the black holes. Thus, the black holes will become 
larger, and this results in a larger amplitude of the grav- 
itational waveform. 



VI. SUMMARY 

In summary, we have constructed from scratch a new 
numerical code based on the BSSN formalism and the 
moving puncture technique. In the code, an FMR/AMR 
algorithm is implemented via the GrACE package, a 
fourth-order spatial finite-differencing scheme and an it- 
erative Crank-Nicholson scheme for time integrations are 
applied in solving the Einstein equations. Some adjust- 
ments of the BSSN formulation for the constraint equa- 
tions in [2^ are also examined in the work. We have com- 
pared the alternative adjustments with the conventional 
ones with a static Schwarzschild black hole and found 
that with both of them the black hole can be evolved 
stably and accurately. We then investigated the viabil- 
ity of several gauge choices through the simulation of a 
single moving black hole with and without spin. In ad- 
dition to obtaining results consistent with those of other 
researchers, we found a new gauge choice with which one 
can also simulate a moving punctured black hole well. 
We next tested our code with the head-on collisions of a 
BBH system in a vacuum and with the perturbation of 
a massless scalar field. The gravitational waveform ob- 
tained in this code from the collision in vacuum is quan- 
titatively the same as that obtained in the work by other 
groups. The purpose of the head-on collision perturbed 
by a scalar field is to understand qualitatively the effect 
of matter on the evolution of binary black holes, as well 
as testing the code further. The result shows that, with a 
specific configuration, the existence of a scalar field could 
delay the merger of binary black holes, as people have ex- 
pected. The strength of the scalar field will significantly 
affect the gravitational waveform. 

The main goal of this work is the construction of a new 
code for the study of numerical relativity. However, re- 
investigating the conventional methods and exploring the 



alternatives to the methods are also emphasized during 
the development of this code. As one can see in Sec. lV Al 
both adjustments for the constraint equations give sta- 
ble and convergent results. However, Fig.[2]also gives the 
impression that the alternate adjustment performs better 
than the conventional adjustment. This simply indicates 
that there is still an issue of the optimal choice of the con- 
straint addition in the Einstein equations. For the gauge 
condition, our study shows that a new gauge choice, i.e.. 
Gauge VII, is able to pass all of our tests. Therefore, this 
gauge could be a possible choice for BBH simulations, al- 
though more investigation is needed. The waveform ob- 
tained from the head-on collisions shows that with the 
code one can simulate the evolution of a BBH system. 
This enables us not only to continue studying numerical 
relativity in more complicated scenarios, like the inspi- 
ral of binary black holes, the recoil problem, and the 
final spin problem, but also to verify the existing results 
and to go even farther. We believe there will never be 
too much double-checking of existing achievements, and 
any further investigations based on them. The result of 
the head-on collision perturbed by a scalar field gives us 
some insight into the possible distortion of a gravitational 
wave, since matter is abundant in most of the gravitating 
systems. 

All of these results verify that the code is reliable and 
ready to be engaged in the study of more realistic, astro- 
physical scenarios and of numerical relativity. We plan 
to use this new code to work on the inspiral of BBH 
systems in the very near future. One black hole can be 
determined totally by only 2 parameters, i.e., mass and 
spin. Thus it is interesting to ask how to determine the 
product black hole of a BBH system from the informa- 
tion of the two initial black holes. The spin expansion 
method in 46] shows us some hint as to how to solve this 
problem, and we plan to study this problem with the new 
code next. 
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APPENDIX A: THE APPROXIMATE 
BOWEN-YORK INITIAL DATA 

In the puncture method, the Hamihonian constraint 
equation of the Bowen-York initial data reduces to 
Eq. ((3T|) . When the hnear momentum P and the spin 
S" of a black hole are small, we can solve the equation 
approximately as 

{dl + dl + dl)u = ^\k^=k,,{i + ^r'. (Ai) 

We chose the coordinates such that P points along the 
z-direction and S lies on the xz coordinate plane, i.e., 
S — {Sx, 0,Sz)- In spherical coordinates, we have 

P'' = P cos 61, 
P^ = -Psm0/r, 
P^ = 0, 

S"' = Sj;C0S(psm9 + SzCosO, (A2) 
— {SxCOSipcosO ~ SzSm0)/r, 
= —SxSViiLp/rsmO. 

Note that trOtp — "^^ sin0 in spherical coordinates. From 
Eqs. (HD and dSOD, k,^ reads 

3P 

Krr = — COS 9, 



Kr 



Kb 



— smO -bx sm(/3, 

Ir 

3sin0 



3P 



(Sx cos i-p cos 9 — Sz sin 9), 



cos 9, 



(A3) 



K, 



3P 



K^^ = — ^ cos 9 sin^ 0. 
Then 

QP^ 5 2 
i^^^-if,, = _[- + -(3cos2^?-l)] 



9 „.,A 1 
185'2 2 1 



^Siq + ^(3cos2 - 1) - sin^ 0cos2^] 



(3cos2 6i - 1)] 



r6 ^3 3" 

+ —PSx sm9 sinip g-S'^S'^ cos sin sin 

(A4) 

By introducing the two new variables, fip = P ■ n/P and 
^is — S ■ n/S, Eq. (jA4p can be rewritten as 

18jS'^ r2 1 ,^ 9 , T 18 , ^ 
+^f3-3(^^^"^)] + 7^(^^^^-"- 

(A5) 



Since Eq. ([3T|) is mainly linear with respect to u, we can 
solve it by dividing K^^ Kij into three parts, i.e., the P^ 
part, the 5^ part, and the P x S part. With a variable 
separation method, we can get the solution of Eq. (pT]) 
[ai,!!!], i.e., Eq. ([32]). 

For the scenario that P — while the spin is large, 
the solution of Eq. (|3T|) can be approximated with the 
spherical symmetric solution [36]. Then Eq. ([3T|) can be 
smeared out as 



In spherical coordinates, this equation reads 



+ -dr)4' 



2r^ 



(A6) 



(A7) 



Following Dain [40| we assume that V' behaves accord- 
ing to a power-law (V' = Ar") and substitute this into 
Eq. (|A7|1 . We find that the approximate solution for a 
highly-spinning black hole is 



(65^) 



2U/8 



(A8) 



However, this approximate solution is valid only for 
rr? I S < r < S/ra, which implies that the larger the 
spin is, the better the approximation is. 



APPENDIX B: WAVE EXTRACTION WITH *4 

In our code we use the Newman-Penrose scalar ^4 for 
extracting the physical information from numerical simu- 
lations. Here we follow Pretorius' work [3] and construct 
the tetrad from the timelike unit normal field as well 
as applying the Gram-Schmidt orthonormalization pro- 
cedure to the three-dimensional vectors 



V = [-y,x,Q\, 

w = [xz.yz, -{x^ + y^)]. 

The tetrad vectors are then given by [ll[ 



n' = ^^, n = ^( u 

V2a V2 a 



(Bl) 



1 



' + (B2) 

771° = 0, m* = —i={v^ + iw^)- 
V2 

From the definition of ^'4 and Gauss-Codazzi-Mainardi 
equations [13, SI] , we have 
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All the quantities in the last equation are physical ones in 
the three-dimensional hypersurface. Let us relate them 
to the conformal ones 



— „40 



e"^{Rijk£ + 'i'y{[iDj]Dk(l} - 4:jk[iDj]De 
160,[i7j][fc0,£] - 8lkiaj]el"'"'(f>,m(f>,n), 



Rijkl 

i?y = R^j -2D,Dj(t>-2%D''Dk(t> 



D^K,k 



(B4) 



3 



Combining these equations together, we are able to 
compute ^'4 on the entire three-hypersurface. When 
we transform to'^ to e*^m^, ^'4 will be transformed to 
e-2«Cvi/4 which results from Eq. (jB3[) . The functions 
with this kind of property are said to have spin weight 
-2. Since the functions with different spin weight are or- 
thogonal, we can calculate the corresponding coefficients 
Airn with respect to the spin weighted spherical harmonic 
function with spin weight -2 [ll!,|4^, i.e., Y^^- The major 
reason for choosing this type of function set is to reduce 
numerical error during the calculations. The coefficients 
of the gravitational radiation can be obtained from the 
integration of 5*4 and each 



*4)=/ / r,-2*4sinM0d(^, (B5) 
Jo Jo 



where £ > 2 and —I < m < £. The spin-weighted spheri- 
cal harmonics can be defined in terms of the Wigner 
d-functions (e.g., [ll|) as 



(B6) 



where 



o-m2{ ) 2^ {i + m-ty.{i + 2-ty. 



t=Ci 



V(£ 2)!(£ + 2)! e . e 

{t-2-m)W. ^2' ^2' ^ ' 

with Ci ~ max(0, m — 2) and C2 — nim{£ + m,£ — 2). 



We have imposed these three "secondary" constraints in 
the calculation of the conformal Ricci curvature, Rij and 
obtain its expression as in Eq. (fTO|) . The explicit expres- 
sion of the scalar curvature R is 



R = -\y'i''%3m + f ^fc + f ^^■'^ (2f,»fe + f,jk ] ■ (C2) 



Here the constraints (|C1[) have been used again. 

The second spatial derivative of the unimodular deter- 
minant constraint gives 



7 lke,ij — li- ku^ J 
and the trace of this "tertiary" constraint is 

7^'7"7.,.fe^ - 4f,,fef(^^)'= = 0. 
Using Eq. (|C4|) we can rewrite Eq. (jC2p to be 



(C3) 
(C4) 
(C5) 



Here the conformal trace curvature is only dependent on 
the first derivatives of the conformal connection and 
of the conformal metric "fij . 

The modified ADM mass volume integral in [2^ 



M 



1 

16^ 



-r^'^f j,fe + (1 - e'^)R 



(C6) 



can be understood better now with the expression (jC5 
Applying Eq. jCH to Eq. jCS]) we have 



M 



1 

16^ 



-ar -e'I'R 



(C7) 



For the cases in which f * is negligible, the ADM mass 
integral expression returns to 



M 



167r 



d'^x 



IQtt p + A,jA'i - -K' 



i'^R 



(C8) 



APPENDIX D: SCALAR FIELD IN A CURVED 
SPACETIME 



APPENDIX C: THE APPLICATION OF THE 
UNIMODULAR DETERMINANT CONSTRAINT 

In the BSSN formulation the determinant of the con- 
formal metric is unimodular, i.e., Eq. (I18|) . The first 
spatial derivative of Eq. gives 

f'l^J,k=0 =^ T'jk=0. (CI) 



The dynamics of a complex scalar field in a curved 
spacetime is described by the following Lagrangian den- 
sity 0, M 

C = -^R+^[V^'^V^<^ + Vm% (Dl) 

where R is the scalar curvature, $ is the scalar field, 
<I> is its complex conjugate, and y(|$p) is the potential 
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depending on Here we consider the case where the 
potential takes the form = m^|<i>p with m being 

the mass parameter. The action with the Lagrangian 
leads to Einstein's equation, 

R^.. - ^gt^uR = SnTf,,, (D2) 

where i?^,^ is the Ricci tensor and T^^ is the stress-energy 
tensor defined by 

T^, = V(^$V,)$ - ^5p.[V,$V^¥ + m2|$|2]. (D3) 

The matter source terms in the BSSN formulation are 
constructed from the stress-energy tensor via the expres- 
sions ilfli] 

p = n^7i.r^", (D4) 
s^J = 7.^7^.^'"^, (D5) 

S^ = -7^pn,^''^ (D6) 

s = 7''sy. (D7) 
The field equation for a massless scalar field is 

V^V^^ = 0, (D8) 



which can be expanded as 

V^V* = -La^(x/^.g'"'5,$) = 0. (D9) 



With the metric given by 

ds^ = -a^dt^ + -fij{dx' + (3'dt){dx-' + (3^dt), (DIO) 

we can re-write the field equation and decompose this 
second-order equation into two first-order-in-time equa- 
tions: 

(dt-p'd^)^ = an, (Dll) 

+ e-^'I'f^a^.^j+aKU. (D12) 

Combining Eqs. (jPlip . (|D12p with Einstein's equation, 
we can study the dynamical interaction between a mass- 
less scalar field and the curved space-time where the 
scalar field lives via the numerical analysis. 
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